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o : 

_ We have examined a metliod of direct extraction of accurate tight-binding parameters from an 

■ ab-initio band-structure calculation. The linear muffin-tin potential method, in its full-potential 
implementation, has been used to provide the hamiltonian and overlap matrix elements in the 
momentum space. These matrix elements are Fourier transformed to real space to produce the 
tight-binding parameters. The feasibility of this method has been tested on the intermetallic alloy 
NiAl, using spd orbitals for each atom. The parameters generated for this alloy have been used as 
input to a real-space calculation of the local density of states using the recursion method. 
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I. INTRODUCTION 



There has been a growing interest in recent years to construct methods of electronic-structure calculation in which 
the required computational time scales linearly with the size of the system, the so-called 0{N) method. [Q The 
proposed approaches abandon the k-space method used in convential band-structure calculations in favor of working 
directly in real space. [^-^ Instead of calculating the energy bands, one instead focuses on local properties, e.g., the 
local density of states. Extensive properties, like the total energy, are obtained by integrating over space instead of 
over the Brillouin zone (hence the linear scaling with the size of the system). This real-space approach is particularly 
relevant to cases in which lack of perfect crystalline symmetry is an essential part of the problem (e.g., in surfaces, 
\ impurities, defects, and amorphous materials) since in those cases Brillouin zone and energy bands are simply non- 
' existent. |^ The tight-binding (TB) method is among the most popular implementations of this real-space strategy. 
O \ The tight-binding approach to electronic structure has a long history which dates back to the classic work of 
Slater and Koster half a century ago. pc|-p^ It has been found especially useful in the modeling of the important 
technological material silicon. jl4-17| The conventional means to obtain the TB parameters is to follow the original 



procedure of Slater and Koster namely to fit the TB energy bands to the ones obtained from accurate self-consistent 
calculations. The parameters for elemental solids obtained in this way have been compiled by Papaconstantopoulos. 
p8[ This database has also been actively enlarged to include parameters for various alloys. p9|-|2^ Recent interest in 
f^F) real-space electronic structure methods has supplied a renewed vigor to this field of TB parametrization. p3|-p7[| 
■ Despite its seemingly simple principle, the fitting process that one uses in obtaining the TB parameters is not 
\ always straightforward in practice. Starting from a guess set of TB parameters one calculates the TB energy bands 
which are then compared with the accurate bands from a self-consistent calculation. One can proceed to minimize the 
merit function by using standard nonlinear optimization procedures, This, in general, is quite a straightforward 
"j^ ■ procedure for simple systems with one atom per unit cell. However, for more complicated multiatom unit cells with 
more intricate "spaghetti" of bands, the number of independent parameters on which the merit function depends 
I I grows rapidly. In the optimization process for this case, the merit function can get trapped easily in a local minimum 
' O ' and the output TB bands have little resemblance to the bands that they are intended to fit. Getting around this 



problem by choosing a good initial set of parameters, unfortunately, still comprises more art than science. Another 
drawback of the TB parameters that are obtained from fitting is that one normally fits only the bands below the 
Fermi energy and a few bands above it. In general the parameters obtained depend on the number of bands that are 
used in the fitting. This feature makes it difficult to attach a physical meaning to the parameters. 

Andersen and coworkers have proposed a method to obtain TB parameters directly from his linear muffin-tin 



?H orbital (LMTO) method. ||29|-|33| This method has also been extended by other groups. |35[| The method outputs the 
. . orthogonal hamiltonian (hopping) parameters between highly-localized (screened) muffin-tin orbitals in the two-center 
approximation. Since each hopping parameter effectively vanishes beyond second-nearest neighbors, TB calculations 
using these parameters can be highly efficient. Moreover, the obtained parameters are also independent of structure 
(transferable), a property which is highly desirable in TB applications. Most of these nice properties, however, depend 
on the atomic sphere approximation (ASA) which is used to derive these results, This approximation has been 
known to work rather well for close-packed structures but is not very accurate for open structures, although this can 
be remedied to some extent by using the so-called empty spheres. The TB-LMTO method assigns a minimal base 
(at most one s, three p, and five d orbitals per atom) to each atom in the crystal. In the LMTO-ASA method, the 
muffin-tin orbitals are set to have zero kinetic energy {k? = 0) in the interstitial region. 



1 



Full-potential (FP) implementations of linear electronic-structure methods, e.g., the FPLMTO method or the FP 
Linear Augmented Plane Wave (FPLAPW) method, are among the most accurate electronic-structure methods. 
p6|-^ It is therefore desirable to have a direct method of extracting practical TB parameters from these methods. 
The basis functions in the FPLAPW method are solutions to Schrodinger equation inside the MT sphere and plane 
waves in the interstitial region. Since the basis functions are not localized, its matrix elements therefore cannot be 
mapped into the TB form. (It should be noted, however, that it is possible to project a FPLAPW calculation into 
localized basis. [^) The basis functions in the FPLMTO method, on the other hand, can be chosen (by setting the 
value of the kinetic energy k^) to have a decaying tail outside its MT sphere. Direct mapping to TB parameters is 
therefore possible for the FPLMTO method. To our knowledge this direct approach has been seriously studied only 
by McMahan and Klepeis who used it to study Si and Si/B phases. p2|-^ 

In this paper we use this direct approach and extract TB (hamiltonian and overlap) parameters for the intermetallic 
alloy NiAl from a FPLMTO method. These parameters are then used as input to a real-space calculation of the local 
density of states using the recursion method. To our knowledge, this is the first test of this direct method to a 
transition-metal alloy. Our choice of NiAl is determined by our previous experience with it and by the fact that 
NiAl has a simple B2 structure in which each atom in the unit cell has the full cubic point group. The latter is of 
practical relevance since the symmetry allows us to reduce the number of independent TB parameters that need to 
be calculated and saved in the database. Since the TB parameters that we extract from the FPLMTO method are to 
be used in a recursion calculation, we do not need to extract the (two-center) Slater-Koster parameters as was done 
by McMahan and Klepeis. j4|] Our motivation in performing this calculation is twofold. First, one hopes to obtain a 
general idea of the level of numerical accuracy that can be expected to result from a combination method like this. 
Second, eventually one would like to use the TB parameters extracted from this method in applications where it is 
too expensive computationally to use direct ab-initio methods, e.g., in studying the effects of a small concentration 
of impurity atoms to static and dynamic properties of alloys. Since the TB method is a non-self-consistent method, 
the plausibility of its result can only be judged by how well it reproduces certain "benchmark" results. In this paper 
the "benchmarks" are the local density of states that can also be obtained from the FPLMTO method. 

Section |l| of this paper discusses the specific implementation of the FPLMTO method that we use for this work. 
The choice of basis is of the utmost importance in a TB calculation. Since the FPLMTO basis functions are well- 



known for being complicated we have chosen to discuss it in some details in Sec. II A . Practical points regarding the 



transformation to real space are listed in Sec. [IB. The recursion calculation and the results that we obtained from 



it are discussed in Sec. [II; and major points of the paper are discussed in a summary section at the end of the paper. 



II. FPLMTO METHOD 

In this paper we have used an implementation of the FPLMTO method that is developed by Wills and Price. 
p6|"H . Several other implementations are available in the literature and, although the main formalism is largely the 
same, each implementation is different in its details. The FPLMTO uses a fixed basis set in calculating its 

hamiltonian and overlap matrix elements; it is therefore similar to the linear combination of atomic orbitals (LCAO) 
method. |^2| However, while in the LCAO method the basis set is determined, and fixed, at the beginning, the basis 
set in the FPLMTO method is determined self-consistently. The FPLMTO method is very similar to the Augmented 
Spherical Wave (ASW) of Wilhams, Kiibler, and Gelatt. Q 

One of the strengths of the FPLMTO method is that it allows unlimited number of basis functions to be assigned 
to each atom in the unit cell. This provides unlimited variational flexibility in finding the ground state configuration. 
The practical requirement of TB, however, requires us to use a minimal basis set. Furthermore, although FPLMTO 
method allows the use of localized (k^ < 0) or unlocalized (k^ > 0) basis, only the localized basis should be 
used for TB purposes. Williams et al. have used ~ —0.2 Ry exclusively in their ASW paper. They have also 
warned that the use of a single fixed n is rather restrictive and cited the work of Gunnarsson et al. [55| that shows 
the variation of eigenenergies with k in a molecular calculation. 

A. Basis Functions 

In the FPLMTO method the crystal is divided into non-overlapping muffin-tin (MT) spheres surrounding each 
atomic sites. The radius of the MT sphere that is centered at r^, the position of the a-th atom in the unit cell, 
is denoted by Sq. On this muffin-tin geometry we construct a muffin-tin potential. This potential is spherically 
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symmetric within each mufHn tin and equal to a constant Vq, the muffin-tin zero, at the interstitial region. The 
potential inside the MT, VMT{f), does not have to join continuously to the interstitial potential Vq at the MT radius. 

We emphasize that in the FPLMTO method the MT potential is used solely to construct the basis functions. The 
final output charge density and potential, which are obtained from a self-consistent calculation, do not necessarily 
have the form of a MT potential. During each iteration in the self-consistency loop, the LMTO basis functions are 
used to solve the Schrodinger (or Dirac) equation in a variational manner similar to the LCAO method. The resulting 
eigenf unctions are then used to calculate the new charge density and potential. The important difference from the 
LCAO method is that here, in addition to the charge density and potential, one also updates the basis functions on 
every step. This update is performed by constructing a new MT potential from the output potential. The constant 
potential Vq is given the value of the average potential at the interstitial, while the spherically-symmetric potential 
i'mt('') is obtained from the angular average of the potential inside the muffin tin. p9| , |5^ -|58[| Note that we average 
the potential rather than the charge density as used in an alternative scheme. |^ 

The final output MT orbitals therefore provide the "best" set of basis functions (within the variational freedom 
imposed by the choice of the LMTO parameters) for performing an LCAO-like calculation on the system under 
consideration. This view of the MT orbitals as basis functions for an LCAO-type calculation has appeared from time 
to time in the literature, under the name of the linear combination of muffin-tin orbitals (LCMTO) method. It has 
been particularly useful for calculations of the electronic structure of systems that do not possess lattice periodicity, 
e.g., surfaces, impurities, and atomic clusters. |]60|-|63|| Unfortunately, the difficulty in performing the necessary three- 
dimensional integrations in this method has made it less popular relative to other methods that utilize orbitals with 
simpler analytical properties like the Slater- or Gaussian-type orbitals. ||56| , |6^ , |6^ 

Each basis function in the FPLMTO method is centered around an atomic site that we will call its host site. With 
respect to this site, all other atomic sites will be referred to as remote sites. The basis function is constructed as a 
continuous patchwork of three different parts: (1) inside the host MT sphere; (2) at the interstitial region; (3) inside 
all remote MT spheres. Thus the basis function centered at atomic site Tq, inside the unit cell with lattice vector 
has the form: 

(l)iair) = (j)ia{ia\r) + (l)ia{I\r) + 0ia(i'a'|r)- (1) 

i' a' 

Here the primed summation means that it should be carried out over all {i'a') ^ (ia). The MTO also carries other 
indices: L = (Im) which are the angular-momentum quantum numbers of the part of the MTO inside the host sphere 
(j)ia{ia\r); the tail parameter n which controls the behavior of the orbital at the interstitial region; and possibly also 
the principal quantum number n which controls the number of radial nodes of (l)ia(ia\r). The last is usually not 
necessary since one normally assigns only one principal quantum number for each (Im). In any case, for clarity we 
will suppress these indices unless their presence is really necessary. 

An MTO is essentially an augmented spherical wave ( ASW) . This means that it is constructed by defining it to 
be a spherical wave (spherical Bessel, Neumann, or Hankel function, i.e. solution of Helmholtz equation in spherical 
coordinates) in the interstitial region: 

0,„(/|r) = Ki{K,na.) ■ i'YUr^a.) ■ e(/|r), (2) 

with a radial part: |4^j5^ 



Here r^Q, = (r — Ria), Ria = (R-i + '''a), fia = \^ia\j and Via 

r is in the interstitial region and otherwise. For < 0, k 
of the spherical Hankel function of the first kind: 



= ^ia/fia- The masking function 6(/|r) is equal to 1 if 
= this spherical wave is usually expressed in terms 



^ [ |k|2(1 + 3z-i +3z-2), 1 = 2, 

where z = j^jr. Since this spherical wave is the envelope of the full MTO we see that, for negative , the MTO has a 
tail that decays exponentially with distance from its host site. The decay rate is controlled directly by the magnitude 
of K. This feature makes it suitable for tight-binding calculations. 



3 



Inside the MT spheres, the envelope spherical wave is augmented (i.e. replaced) with a solution of the Schrodinger 
equation for the MT potential. The augmentation process is similar to the augmentation of the plane wave in the 
Augmented Plane Wave (APW) method. 1^,^ The difference, of course, is that here we are augmenting a spherical 
wave instead of a plane wave. Note that the basis functions in the APW method are not suitable for tight binding 
since their plane-wave envelopes are not localized in space. 

As we approach the host MT sphere from the outside, the basis function is determined by Eq.(|^). Note that, relative 
to the origin at the host site, this function is separable, i.e. it is a product of a radial and an angular part. Inside the 
host MT sphere, this function is augmented with a solution of the Schrodinger equation for the MT potential. Since 
this potential is spherically symmetric inside the MT, the basis function inside the host MT is also separable: 



Ua r I 



[Aipiai (e 

nl ■} ^ia ial {^nl 1 '^ia 



(5) 



where the MT masking function Q{ia\r) is equal to 1 if r is inside the MT sphere (ia) and otherwise. The energy 
dependence of the radial solution has been approximated with a linear combination of the radial solution tpiai and 
its energy derivative (piai (= dipiai/de) calculated at a fixed chosen energy e„;. |Q The coefficients are obtained by 
requiring the function to match continuously and smoothly with the interstitial function at the MT radius: 





'A' 




KiiK,Sia) 




B 




if J (k, Sia ) 



(6) 



Here a prime denotes differentiation with respect to the radial coordinate: ip' = dtp/dr. 

The augmentation process at the remote spheres proceeds rather differently from the one for the host sphere. As 
we approach a remote MT sphere from the interstitial region, the basis function is again determined by Eq.(^. The 
important difference from the host-sphere case is that relative to the origin at the remote site the interstitial function 
is not separable. We can, however, express the spherical Hankel function at the interstitial region as a sum of separable 
hmctions in the neighbourhood of the remote site: |65| - |69[| 



(7) 



with the spherical Bessel function as the radial part in each term: 

sinh z, 



Ji(^K,r) = K ''ji{Kr) = - X ^ \k\ -^(coshz — z -^sinhz), 

^|-2[(3z- - 



I = 0, 

l)sinhz — 3z ^coshz], 1 — 2. 



(8) 



The addition theorem for the Bessel functions, Eq.(0), is a consequence of the fact that both sides of the identity 
are solution of the translationally- invariant Helmholtz equation. The spherical Hankel function, Ki{K,r), is regular 
everywhere except at the origin, while the spherical Bessel function, J/(k, r), is regular everywhere except at infinity. 
The domain of validity of Eq.(Q) is just the domain where both sides of the identity are regular. 

The radial part of each term in the one-center expansion, Eq.(^, is centered at the remote site R^'q'. Thus the 
{(fijif) augmentation of the interstitial function, Eq.(|2|), at this remote site can be performed by augmenting each 
Jii^K, Vila') with the radial part of the solution of the Schrodinger equation for the spherically-symmetric MT potential 
centered at Rj'a'- Thus we replace: 



(9) 



The coefficients are again obtained by requiring continuity in the value of the function and its first derivative with 
those of the interstitial envelope function at the MT radius: 





'ci 












Jj^ (/^, Si' (x' ) 



The part of basis function inside a remote MT sphere is thus given by: 

Kim 



,(z'a'|r) = e(i'a'|r) • ^ [Q'lpi'a'i'in'a') + Di^tpi'a'i'in'a')] • ^l' (r^'a') ' Ri'a' - R^a)- 



(10) 



(11) 
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The summation on the r.h.s. of the equation ideally runs to infinity. Practically, a converged total energy is achieved 
using / < ~ 6 — 8 in most cases. Note that, inside the remote sphere, an angular-momentum expansion is 
necessary since the resulting function needs to match the envelope interstitial function at its MT sphere. Relative to 
this remote site, the envelope function contains multiple harmonics since it is centered at a different site. 

The MT orbital, constructed from Eqs.(|^), and (11), is continuous and smooth everywhere. Note, however. 



that it is not normalized and it does not have a separable form: 

(j)ia{r) ^ u{ria)Yi,n{ria), (12) 

in other words, the MTO does not have good angular-momentum quantum numbers. Although the parts of the 
MTO inside the host MT and at the interstitials are separable, the presence of the remote sites breaks the perfect 
spherical symmetry around the host site, and therefore also removes the separability. The MTO, however, does obey 
exactly the point symmetry of its host site. Moreover, since the non-separable contributions come from the remote 
sites, where the amplitude of the exponentially-decaying envelope function is small, to a good approximation the 
MTO does transform as if it has the quantum numbers (Zm) of the part of the orbital inside the host sphere. The 
fact that the MTO is not fully separable should be noted especially if one intends to map the FPLMTO result into 
a TB form using the Slater-Koster parametrization method. In the SK-LCAO parametrization method the basis 
functions are assumed to be atomic orbitals; these orbitals are separable since they are derived from a single-atom 
spherically-symmetric potential. [|lO|- 



B. Fourier Transform 



In applications of the FPLMTO method to ideal crystalline systems, instead of using the basis functions in real 
space, one works with the Bloch functions: 



,(k,r) = ^ exp(ik • Ri) <j>ia{r). (13) 



The energy bands are obtained by solving the eigenvalue equation: 

det\H^p{k)~ESc.fj{k)\^0. (14) 

The rank of the matrices is equal to the total number of the orbitals used for all the atoms within the unit cell. The 
hamiltonian matrix is given by integration over an arbitrarily chosen unit cell R„: 

H^k) = (V'a(k,r)|iJ|V'Mk,r)) = l^e^k.(R,-R,)^ jj |^,0(/„)) -|- ^(0,„(n7)| H |</.,M"7))], (15) 

ij n 7 

and similar expression is obtained for the overlap matrix by replacing the hamiltonian operator H with identity. The 
vectors have been defined by: 

{r\(j)ia{nj)) = (l)ia{nj\r), (16) 

which is the part of the basis function, centered at the host site ia, inside the MT sphere 717 or the interstitial region 
In-, which is defined to be the interstitial region within the unit cell R„. Cross terms are absent from the square 
bracket in Eq.(^5|) since the overlap of their masking functions is zero. 

It should be mentioned that in the actual FPLMTO formalism, the matrix elements are not computed using a 
multicenter expansion (over all lattice sites) as shown in Eq.(p^. Rather, using the addition theorem for the tail 
of the MT orbital, Eq.(^, the multicenter expansion is first transformed into a one-center expansion (over angular 
momenta at a single site), and this sum is the one that is actually computed. jS^,^^ This, however, is merely a 
computational contrivance, albeit a pivotal one in the FPLMTO method, and should not mask the true structure of 
the basis function in k-space as a simple Bloch sum of the MT orbitals. 

In the FPLMTO method of Wills et al. [Q , the fact that each term of the matrix elements is to be calculated only 
over a certain region, either inside the MT sphere or the interstitial region as seen in Eq.(p^, is used to facilitate its 
computation. Instead of using the true MTO, which is complicated, one uses a pseudo basis orbital. |Q This orbital 
is equal to the true MTO inside the relevant integration region but outside the region it is replaced with a smooth 
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function which is chosen to facilitate the computation. The integration can then be performed efficiently by working 
in the momentum space. 

Using the translational properties of the Bloch functions, the matrix element can be written as: 

H^pik) = ^exp(ik • R,) iJa/3(Rj), (17) 

j 

with the matrix elements in real space: 

H^^CR,) = [iMln)\ H \c^jp{I^)) + Y,{Mn-i)\ H \(t>jp{n-f))] . (18) 



To calculate the inverse transform of Eq.(17) we use the fact that the real-space matrix element HapCR) decays 
exponentially with the distance |R|. It is thus reasonable to ignore matrix elements between orbitals separated by a 
distance greater than a certain cutoff radius Re- In practice this can be performed by working with a finite crystal 
which contains unit cells along the a;-direction such that N^a^ > 2i?c, with being the lattice constant along 
the a;-direction, and similar ranges for the y and z directions. The real-space matrix elements can then be obtained 
from a discrete Fourier transform: 

H^piRj) = j^^expi^tk ■ R,) Hap{k). (19) 

k 

Here N = N^NyN^ and ki = imi/ [Nitti) with —Ni < Ui < Ni for i e {x, y, z}. 

Some practical remarks should be mentioned regarding this calculation of the inverse Fourier transform. First, the 
output matrix elements from the FPLMTO are calculated using the spherical harmonic i''Yim{v) as the angular part of 
the basis function inside its host sphere. Usually it is easier to work with real basis functions and this can be achieved 
by making linear combinations of the spherical harmonics. fn\ Second, the basis functions are not normalized. The 
normalized matrix elements are obtained by using the overlap matrix elements: 

<^(R) = 5„a(0)-5 . i/„^(R) . SppiO)-^. (20) 

Third, the actual displacement vector connecting the centers of the two orbitals in HapCR) is not R but is instead 
(R + Ti3 ~ Ta). Instead of the original hamiltonian matrix in Eg. (p^ , it is preferable to work with a modified matrix 
that is related to the original by a unitary transformation: 

i/^'^(k) = e^^^^-^-^H^pik) = J2 e'''-^^^ i?a/3(Rj)- (21) 

j 

This positions the origin at the center of orbital (Oa) and allows the matrix elements to obey the point symmetry of 
the corresponding site. It may therefore be used to reduce the number of matrix elements that need to be stored. 
Note that the eigenvalues are unaffected by this unitary transformation. 



C. Lowdin Orthogonalization 

The real-space hamiltonian matrix elements in Eq. (^9|) (together with the overlap matrix elements Sap which can be 
obtained similarly) can be used directly in a nonorthogonal recursion calculation. [fz^-fT^ Unfortunately, this requires 
us to obtain the inverse of the overlap matrix in real space, which is a non-trivial computational task. In this paper we 
instead use an orthogonal hamiltonian which is obtained by Lowdin orthogonalization of the original hamiltonian and 
overlap matrices. |0 An alternative scheme is to use the chemical pseudopotential approach and work with a non- 
hermitian matrix S~^H. [ |75| , [76t The advantage of Lowdin orthogonalization is that it is a symmetry transformation: 
angular symmetry is preserved, i.e., orbitals with angular momentum / are not mixed by the transformation with 
orbitals with V ^ I. This means that the real-space matrix elements can still be parametrized by Slater-Koster 
parametrization, e.g., using the procedure proposed by McMahan and Klcpeis. 

Since we are working with a perfect crystal, the Lowdin orthogonalization process can be applied to the momentum- 
space matrix elements. One first diagonalize the overlap matrix: 

D = U^ -S-U, (22) 
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FIG. 1. Energy bands for NiAl (only the 9 lowest bands are shown). The lines are the output eigenvalues from FPLMTO 
while the dots are calculated using the TB parameters. The two sets coincide with each other since the TB parameters are 
obtained by direct extraction from the FPLMTO. The horizontal line is drawn at the Fermi energy, Ef = 1.0417 Ry. 



where U is the matrix containing the eigenvectors of S as its columns, and D is a diagonal matrix containing the 
eigenvalues of S. Using these matrices, we define a new matrix: 

A = U-D-i-Ul (23) 

Since the overlap matrix is a positive-definite hermitian matrix, its eigenvalues are all positive real numbers and the 
matrix D^^ is well defined. In practice, however, small negative eigenvalues of S may sometimes appear due to the 
finite machine precision used in the computation. This problem is especially relevant if the diagonalization of the 
overlap matrix is performed in real space and when one uses an overlap matrix which is necessarily approximate due 
to the finite cutoff imposed on the range of the matrix elements or due to the Slater-Koster two-center approximation 
used in parametrizing them. One commonly-used fix to this problem is to add a diagonal matrix with small elements to 
make the overlap matrix positive definite. Jt^ The orthogonal Lowdin hamiltonian matrix is obtained by sandwiching 
the original hamiltonian matrix with A matrices: 

H^^'^^A-H-A. (24) 

The nonorthogonal eigenvalue problem, Eq.(p^), is then transformed to an equivalent orthogonal one: 

det|i/(J(k)-£;(5„;3(k)| =0. (25) 

The orthogonal hopping parameters are obtained by substituting H^^J (k) in Eq. ( [T9| ) . These are the parameters that 
we use in the recursion calculation which is described in the rest of this paper. 



III. RECURSION CALCULATION 



The recursion method, initially introduced to electronic-structure calculation by Haydock, is one of the most efficient 
and stable ways of calculating the electronic Green's function from the TB hamiltonian describing the system: W9 



Guv{E)^{u\{E-H)-^\v), (26) 

where \u) and \v) are two arbitrary vectors in the Hilbert space for the problem. The local density of states (LDOS) 
is a projection of the total density of states to a local orbital: 

Pa(£;)=^|(a|*„>P<5(i?-i?„), (27) 

n 

where jv^n) is an eigenvector with eigenvalue En- This can be obtained from the Green's function: 

Pa [E)^-- Im [G„„ {E + it)], (28) 
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FIG. 2. Recursion parameters for the LDOS of Ni-d orbitals. Filled circles (•) are the a„ parameters while open circles (o) 
are b„ parameters. 
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FIG. 3. Recursion parameters for the LDOS of Al-p orbitals. Filled circles (•) are the a„ parameters while open circles (o) 
are bn parameters. 



with e — > 0"*". The token Im[- • •] denotes the operation of taking the imaginary part of its argument. The recursion 
method outputs the Green's function as a continued fraction: 

GiE) = 1 ^ . (29) 

{E - ao) ^ 



(E-ai) 

Note that there are two conventions used in the literature regarding the indexing of the a-paramctcrs in the above 
continued fraction. One convention is to assign ao as the first a-parameter, [ |7^ , [73[ while the other use ai as the 
first parameter. Here we follow the first convention. The calculation of the LDOS is therefore reduced to the 
computation of the strings of recursion coefficients {a„} and {bn}- The procedure for performing this has been 
described extensively in the literature and will not be repeated here. | |7^ , p0[ 

We have used the recursion method to calculate the LDOS for NiAl crystal. This alloy crystallizes in B2 structure 
with two atoms per (cubic) unit cell. One component (Ni or Al) sits at (0, 0, 0), while the other is at the body-center 
of the cube, ^, ^). We obtained the real-space orthogonal hopping parameters by direct extraction and Lowdin 
transformation as described in the preceding section. The parameters were extracted from a single-K FPLMTO 
calculation using 9 spd orbitals and = —0.04 Ry for each Ni or Al atom. This rather small value of k produces 
basis functions with relatively long range (the envelope of the basis function decays roughly as e"'"''"). Using a larger 
negative value of k, however, produces a greater discrepancy between the output bands from this single-K calculation 
and the corresponding accurate bands from a multiple-K calculation. The chosen value of k that we used is obtained by 
compromising the need to have a short-ranged TB basis functions with the requirement to have a precision comparable 
to accurate multiple-/* result. 

All of the extracted parameters are saved in a database for use in the recursion calculation (crystal symmetry was 
used to reduce the size of the database). They were calculated for separation distance of up to 8 lattice constants. The 
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FIG. 4. The local density of states (LDOS) for Ni and Al states in NiAl as calculated using the tight-binding recursion method 
(dark full line) and from Brillouin-zone integration of FPLMTO eigenvalues with atom and angular-momentum projection (grey 
line with plus (+) signs). (a,b,c) The left column shows s,p, d LDOS's for Ni; (d,e,f) the right column displays the a,p, d LDOS's 
for Al. The Fermi energy is at Ef = 1.0417 Ry. 



9 



obtained TB parameters fall roughly one order of magnitude for each increase in separation distance of one lattice 
constant. The recursion calculation was performed using a cubic cluster of (16)'^ unit cells and the LDOS's were 
computed for Ni and Al orbitals within the unit cell at the center of the cluster. Since the TB parameters connecting 
the orbitals at the central unit cell to the ones at the surface of the cluster are practically zero, the calculated LDOS 
should not depend on the size of the cluster if we increase it further. 

Figs. H and ^ show the calculated recursion parameters for Ni-d and Al-p states, respectively. Here the d-state is 
defined by: 

1 2 

= E M2,m), (30) 

with a corresponding definition for the p-state. It is seen that the parameters tend to settle around certain constant 
values. Indeed it is known that as n — > oo: 

a„ ^ {Et + Et)/2, 6„ ^ W/4, (31) 

where Et and Eb are the values of energy at the top and bottom of the spectrum, respectively, and W = {Et — Eh) 
is the total bandwidth. Furthermore, asymptotically the parameters should oscillate around their limit values with a 
decaying amplitude as n ^ oo. The frequency and the rate of decay are determined by specific features (the Van 
Hove singularities) of the spectrum. |Q 

The recursion parameters shown in Figs. |^ and ^ are seen to roughly follow this decayed oscillation prediction. 
However, for large values of n (greater than ^ 40), the amplitudes of the deviation tend to increase rather than 
decrease. We interpret this as the limit at which the machine precision noise start to interfere with the calculation. 
This interpretation is supported by our experience that incorporating the recursion parameters with n > in the 
reconstruction of LDOS in general does not help improve the agreement, and indeed it tends to reduce the agreement, 
with the accurate LDOS obtained from the FPLMTO method. Note that the main workhorse in the recursion method 
involves a matrix-vector multiplication with the dimension of the vector being the total number of orbitals in the 
cluster. In our case this dimension is equal to 73728 (= 16^ unit cells x 2 atoms/cell x 9 orbitals/atom). Compounded 
with this large dimension is the fact that the elements of the matrix and the vector in general have widely different 
orders of magnitude which makes the calculation rather prone to rounding errors. 

Figs. ^ display the main results in this paper: the s,p,d LDOS of Ni and Al in NiAl as calculated using the TB 
recursion method. These LDOS's are compared with the corresponding accurate spectra obtained from Brillouin zone 
integration using atom and angular-momentum projection of the eigenvalues obtained directly from the FPLMTO 
method. The agreement in general is good although the TB recursion method is unable to precisely reproduce the 
sharp peaks in the FPLMTO spectra. Not very surprisingly, the best fit is obtained for the tightly-bound Ni-d state. 

The limited number (n„j ~ 40) of reliable recursion parameters that we can use here is generally not sufficient for 
acceptable reproduction of the FPLMTO LDOS. Furthermore, construction of LDOS by directly taking the imaginary 
part of the Green's function, Eq.(p9|), normally introduces an undesirable feature in the output LDOS in the form 
of spurious rapid oscillations near the edges of the spectrum, as can still be seen in, e.g., Fig. ||a. The simplest 
extrapolation method for the recursion parameters is the square-root terminator method: ||79| ] one simply replaces all 
an and 6„ for n > rim with constants aoo, and hoo- The tail of the continued fraction for the Green's function can 
then be summed exactly as a square-root expression. The best values for the terminating coefficients Ooo and boo can 
be obtained using the Beer-Pettifor method. [Q. Other ways to construct a terminator for the continued fraction 
have been suggested, but for our calculation we have chosen to use the Beer-Pettifor method. Finally, further 
smoothing of the spectra shown in Fig. ^ have been performed using a Chebyshev-polynomial method. Js^ 



IV. SUMMARY 



In this paper, we have described a method for directly extracting real-space TB parameters from a FPLMTO 
method. The basis functions used in the FPLMTO calculation and their construction are described in considerable 
details. Special emphasis has been placed on the fact that these basis functions do not have a well-defined angular- 
momentum quantum number, although to a good approximation they do. We believe this fact should be kept in mind 
in transcribing FPLMTO matrix elements to the Slater-Koster two-center form, since there it is implicitly assumed 
that the basis functions do have a good angular-momentum quantum number. 

This direct extraction method has been applied to an intermetallic alloy NiAl. The TB parameters extracted have 
been used as input to a real-space calculation of local densities of states using the recursion method. We believe 
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this is the first apphcation of the direct extraction method to a real-space calculation of the electronic structure of 
an intermetallic alloy. The good agreement between the LDOS obtained using the TB extracted from the FPLMTO 
method and the accurate LDOS obtained directly from the FPLMTO method shows the feasibility of using this 
method to study more complicated cases, e.g., alloys with small concentration of impiu'ity atoms. 

This work was supported by AF-OSR Grant No. F49620-99- 1-0274. DD would like to thank Dr. M.J. Mehl of 
the Naval Research Laboratory, Washington D.C., for useful discussions and insightful inputs regarding the NRL 
Tight-Binding code. 
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